Time complexity and convergence analysis of domain theoretic Picard method

We present an implementation of the domain-theoretic Picard method for solving initial value problems (IVPs) introduced by Edalat and Pattinson [1]. Compared to Edalat and Pattinson's implementation, our algorithm uses a more efficient arithmetic based on an arbitrary precision floating-point library. Despite the additional overestimations due to floating-point rounding, we obtain a similar bound on the convergence rate of the produced approximations. Moreover, our convergence analysis is detailed enough to allow a static optimisation in the growth of the precision used in successive Picard iterations. Such optimisation greatly improves the efficiency of the solving process. Although a similar optimisation could be performed dynamically without our analysis, a static one gives us a significant advantage: we are able to predict the time it will take the solver to obtain an approximation of a certain (arbitrarily high) quality

Fast linear algebra is stable

In an earlier paper, we showed that a large class of fast recursive matrix multiplication algorithms is stable in a normwise sense, and that in fact if multiplication of $n$-by-$n$ matrices can be done by any algorithm in $O(n^{\omega + \eta})$ operations for any $\eta > 0$, then it can be done stably in $O(n^{\omega + \eta})$ operations for any $\eta > 0$. Here we extend this result to show that essentially all standard linear algebra operations, including LU decomposition, QR decomposition, linear equation solving, matrix inversion, solving least squares problems, (generalized) eigenvalue problems and the singular value decomposition can also be done stably (in a normwise sense) in $O(n^{\omega + \eta})$ operations.Comment: 26 pages; final version; to appear in Numerische Mathemati

New, Highly Accurate Propagator for the Linear and Nonlinear Schr\"odinger Equation

A propagation method for the time dependent Schr\"odinger equation was studied leading to a general scheme of solving ode type equations. Standard space discretization of time-dependent pde's usually results in system of ode's of the form u_t -Gu = s where G is a operator (matrix) and u is a time-dependent solution vector. Highly accurate methods, based on polynomial approximation of a modified exponential evolution operator, had been developed already for this type of problems where G is a linear, time independent matrix and s is a constant vector. In this paper we will describe a new algorithm for the more general case where s is a time-dependent r.h.s vector. An iterative version of the new algorithm can be applied to the general case where G depends on t or u. Numerical results for Schr\"odinger equation with time-dependent potential and to non-linear Schr\"odinger equation will be presented.Comment: 14 page

Robust high-dimensional precision matrix estimation

The dependency structure of multivariate data can be analyzed using the covariance matrix $\Sigma$. In many fields the precision matrix $\Sigma^{-1}$ is even more informative. As the sample covariance estimator is singular in high-dimensions, it cannot be used to obtain a precision matrix estimator. A popular high-dimensional estimator is the graphical lasso, but it lacks robustness. We consider the high-dimensional independent contamination model. Here, even a small percentage of contaminated cells in the data matrix may lead to a high percentage of contaminated rows. Downweighting entire observations, which is done by traditional robust procedures, would then results in a loss of information. In this paper, we formally prove that replacing the sample covariance matrix in the graphical lasso with an elementwise robust covariance matrix leads to an elementwise robust, sparse precision matrix estimator computable in high-dimensions. Examples of such elementwise robust covariance estimators are given. The final precision matrix estimator is positive definite, has a high breakdown point under elementwise contamination and can be computed fast

Analytical study of non-linear transport across a semiconductor-metal junction

In this paper we study analytically a one-dimensional model for a semiconductor-metal junction. We study the formation of Tamm states and how they evolve when the semi-infinite semiconductor and metal are coupled together. The non-linear current, as a function of the bias voltage, is studied using the non-equilibrium Green's function method and the density matrix of the interface is given. The electronic occupation of the sites defining the interface has strong non-linearities as function of the bias voltage due to strong resonances present in the Green's functions of the junction sites. The surface Green's function is computed analytically by solving a quadratic matrix equation, which does not require adding a small imaginary constant to the energy. The wave function for the surface states is given

User-friendly tail bounds for sums of random matrices

This paper presents new probability inequalities for sums of independent, random, self-adjoint matrices. These results place simple and easily verifiable hypotheses on the summands, and they deliver strong conclusions about the large-deviation behavior of the maximum eigenvalue of the sum. Tail bounds for the norm of a sum of random rectangular matrices follow as an immediate corollary. The proof techniques also yield some information about matrix-valued martingales. In other words, this paper provides noncommutative generalizations of the classical bounds associated with the names Azuma, Bennett, Bernstein, Chernoff, Hoeffding, and McDiarmid. The matrix inequalities promise the same diversity of application, ease of use, and strength of conclusion that have made the scalar inequalities so valuable.Comment: Current paper is the version of record. The material on Freedman's inequality has been moved to a separate note; other martingale bounds are described in Caltech ACM Report 2011-0

Accuracy and Stability of Computing High-Order Derivatives of Analytic Functions by Cauchy Integrals

High-order derivatives of analytic functions are expressible as Cauchy integrals over circular contours, which can very effectively be approximated, e.g., by trapezoidal sums. Whereas analytically each radius r up to the radius of convergence is equal, numerical stability strongly depends on r. We give a comprehensive study of this effect; in particular we show that there is a unique radius that minimizes the loss of accuracy caused by round-off errors. For large classes of functions, though not for all, this radius actually gives about full accuracy; a remarkable fact that we explain by the theory of Hardy spaces, by the Wiman-Valiron and Levin-Pfluger theory of entire functions, and by the saddle-point method of asymptotic analysis. Many examples and non-trivial applications are discussed in detail.Comment: Version 4 has some references and a discussion of other quadrature rules added; 57 pages, 7 figures, 6 tables; to appear in Found. Comput. Mat